# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
# install.packages(c("tidyverse"))

library(tidyverse)
# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")

df_amtrak <-
  mturk_opeds %>%
  transmute(Z = as.numeric(Z == "amtrak"),
            Y = dv_amtrak_agree_w1,
            oped = "\"The Amtrak crash: Is more spending the answer?\" \n Randal O'Toole in Newsweek")
df_climate <-
  mturk_opeds %>%
  transmute(Z = as.numeric(Z == "climate"),
            Y = dv_climate_agree_w1,
            oped = "\"The political assault on climate skeptics\" \n Richard S. Lindzen in the Wall Street Journal")
df_paul <-
  mturk_opeds %>%
  transmute(Z = as.numeric(Z == "paul"),
            Y = dv_flat_agree_w1,
            oped = "\"Blow up the tax code and start over\" \n Rand Paul in the Wall Street Journal")
df_veterans <-
  mturk_opeds %>%
  transmute(Z = as.numeric(Z == "veterans"),
            Y = dv_vets_agree_w1,
            oped = "\"The other veterans scandal\" \n Michael F. Cannon and Christopher Preble \n in The New York Times")
df_wallstreet <-
  mturk_opeds %>%
  transmute(Z = as.numeric(Z == "wallstreet"),
            Y = dv_wall_agree_w1,
            oped = "\"Wall Street offers very real benefits\" \n Thaya Knight in USA Today")

gg_df <- bind_rows(df_amtrak,
                   df_climate,
                   df_paul,
                   df_veterans,
                   df_wallstreet) %>%
  mutate(treatment = if_else(Z == 1, "Treatment", "Control")) %>%
  filter(!is.na(Y)) %>%
  group_by(oped, treatment) %>%
  summarize(
    group_mean = mean(Y) * 100,
    group_n = n(),
    group_se = sd(Y) / sqrt(group_n) * 100,
    ci_lower = group_mean - 1.96 * group_se,
    ci_upper = group_mean + 1.96 * group_se
  )


gg_df <-
  within(gg_df, {
    label = c(rep(NA, 8), "Did not read op-ed", "Read op-ed")
    hjust = c(rep(NA, 8), "outward", "outward")
  })

pro_con_colors <- c("#C67800", "#205C8A")

g <-
  ggplot(gg_df,
         aes(
           x = oped,
           y = group_mean,
           group = treatment,
           fill = treatment,
           color = treatment
         )) +
  #geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), position = position_dodge(width = .75), width = 0) +
  geom_bar(stat = "identity",
           position = position_dodge(width = .75),
           alpha = 0.5) +
  geom_text(
    aes(
      y = (group_mean + 2),
      label = label,
      hjust = hjust
    ),
    color = "black",
    position = position_dodge(width = .75)
  ) +
  coord_flip() +
  theme_bw() +
  scale_fill_manual(values = pro_con_colors) +
  scale_color_manual(values = pro_con_colors) +
  ylab("Percentage Agreeing with Op-Ed Author (dichotomized agreement scale)") +
  ylim(0, 100) +
  ggtitle("Percent agreeing with op-ed") +
  labs(caption = "Source: Coppock, Ekins, and Kirby (2018), MTurk sample") +
  theme(
    legend.position = "none",
    plot.caption = element_text(hjust = 0),
    axis.title = element_blank()
  )

g
